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Abstract.  We  study  the  simulation  of  specular  reflection  in  a  level  set  method  imple¬ 
mentation  for  wavefront  propagation  in  high  frequency  acoustics  using  WENO  spatial 
operators.  To  implement  WENO  efficiently  and  maintain  convergence  rate,  a  rectangu¬ 
lar  grid  is  used  over  the  physical  space.  When  the  physical  domain  does  not  conform  to 
the  rectangular  grid,  appropriate  boundary  conditions  to  represent  reflection  must  be 
derived  to  apply  at  grid  locations  that  are  not  coincident  with  the  reflecting  boundary. 
A  related  problem  is  the  extraction  of  the  normal  vectors  to  the  boundary,  which  are 
required  to  formulate  the  reflection  condition.  A  separate  level  set  method  is  applied 
to  pre-compute  the  boundary  normals  which  are  then  stored  for  use  in  the  wavefront 
method.  Two  approaches  to  handling  the  reflection  boundary  condition  are  proposed 
and  studied:  one  uses  an  approximation  to  the  boundary  location,  and  the  other  uses 
a  local  reflection  principle.  The  second  method  is  shown  to  produce  superior  results. 

AMS  subject  classifications:  65M06, 65M22, 65M25 
PACS:  43.30.Gv,  43.30.Zk 

Key  words:  Boundary  conditions,  reflection,  level  set  method,  wavefront  methods,  high-frequency 
acoustic  propagation,  WENO. 


1  Introduction 

1.1  Motivation 

High  frequency  wave  propagation  models  are  of  vital  importance  to  understanding  the 
propagation  of  sound  in  shallow  water  environments.  Applications  such  as  acoustic  to¬ 
mography  [1]  and  underwater  communications  [2]  rely  on  ray  tracing  for  system  design 
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and  performance  prediction.  The  frequencies  useful  for  such  applications  in  shallow 
water  are  at  least  on  the  order  of  1kHz  and  can  range  into  the  MHz  regime,  rendering 
full  wave  models  impractical  in  simulations.  On  the  other  hand,  ray  tracing  has  many 
known  limitations.  In  particular,  the  Lagrangian  nature  of  the  ray  trace  model  leads  to 
difficulty  resolving  wave  arrivals  at  a  fixed  location  in  space.  This  is  especially  true  in 
range-dependent,  shallow  water  environments.  The  goal  of  this  work  is  to  devise  stable, 
accurate,  numerical  boundary  conditions  for  reflection  in  a  level  set  method  implemen¬ 
tation  that  solves  the  high  frequency  equation  for  the  acoustic  phase,  thus  tracing  entire 
wavefronts  on  a  fixed  grid  in  space.  A  practical  wavefront  model  for  shallow  water 
acoustics  would  allow  for  improved  accuracy  in  such  harsh  environments,  where  multi- 
path  propagation  is  evident  due  to  multiple  reflections  off  of  the  sea  surface  and  bottom. 
Such  a  model  would  also  be  more  robust  to  perturbations  in  material  properties.  For 
example,  Godin  [3]  studied  the  effect  of  small  perturbations  in  the  sound  speed  on  rays, 
showing  that  the  resulting  displacement  of  rays  along  the  wavefront  are  much  smaller 
than  displacements  across  the  wavefront.  This  type  of  relative  stability  of  wavefronts  is 
critical  to  practical  simulations  due  to  the  uncertainty  present  in  sound  speed  profiles 
derived  from  at-sea  data  collection. 

This  study  is  based  on  a  method  laid  out  in  [4],  introducing  an  implementation  of  the 
level  set  method  intended  for  a  two-dimensional  shallow  water  acoustics  application. 
The  algorithm  is  based  on  the  level  set  method  for  geometric  optics  [5].  The  level  set 
equations  are  solved  using  a  fifth  order  Weighted  Essentially  Non-Oscillatory  (WENO) 
method  for  the  spatial  operator  [6]  combined  with  third  order  Total  Variation  Diminish¬ 
ing  Runge-Kutta  (TVDRK)  time  integration  [7]  for  stability.  This  choice  is  motivated  by 
the  presence  of  a  discontinuity  in  the  phase  space  at  a  reflecting  boundary  [4].  While  the 
discontinuity  is  a  consequence  of  the  rectangular  scatterers  presented  in  that  work,  this  is 
not  the  case  in  general,  but  even  in  the  benign  case  of  a  flat  reflecting  boundary,  a  sharp 
cusp  forms  at  the  boundary  in  the  reduced  phase  space.  The  examples  presented  in  [4] 
are  limited  by  the  types  of  domain  geometry  that  could  be  handled  under  the  reflection 
boundary  condition  to  either  grid-conforming  rectangular  domains,  or  a  domain  with  a 
linear  boundary.  For  the  linear  boundary  case,  the  boundary  did  not  conform  to  the  grid 
and  an  approximation  was  applied  in  order  to  provide  the  required  boundary  conditions 
for  the  numerical  solver.  This  approach  is  discussed  in  greater  detail  in  Section  2.3.1.  The 
goal  of  the  present  work  is  to  extend  the  method  to  accommodate  more  realistic  domain 
geometry,  and  improve  the  accuracy  of  the  reflection  boundary  condition.  To  achieve  this, 
one  could  modify  the  rectangular  grid  to  conform  to  the  geometry,  as  in  [8],  however  in 
the  general  case,  this  affects  the  convergence  of  WENO.  Refining  the  grid  to  match  the 
geometry  can  result  in  strong  restrictions  on  the  CFL  condition.  Another  approach  would 
be  to  use  finite  elements  or  finite  volume  methods  on  an  unstructured  grid,  but  this  com¬ 
plicates  the  implementation  of  WENO,  especially  in  high  dimensions  (for  a  dimension  d 
physical  space,  the  reduced  phase  space  has  dimension  2d— 1).  Alternately,  the  approach 
taken  in  this  work  is  to  embed  the  reflector  in  the  uniform,  rectangular  grid  and  derive 
appropriate  boundary  conditions  to  apply  at  grid  points  adjacent  to  the  reflector. 
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An  additional  consideration  is  the  knowledge  of  the  boundary  normals  which  is  es¬ 
sential  to  computing  the  reflection  condition.  For  general  domain  geometry,  the  normals 
need  to  be  computed.  We  accomplish  this  via  a  co-dimension  one  level  set  method  that 
computes  the  signed  distance  function  to  the  boundary  [8].  Given  the  signed  distance 
function,  the  normal  vector  can  be  extracted  by  taking  its  gradient,  using  central  differ¬ 
encing,  for  example.  The  procedure  only  needs  to  be  run  once  before  the  main  time  loop, 
and  the  resulting  normal  angles  may  be  stored  for  use  within  the  main  program. 

1.2  Related  work 

In  [4],  a  method  was  described  for  implementing  the  level  set  method  introduced  in 
[5]  for  the  purpose  of  modeling  high  frequency  underwater  acoustic  propagation.  The 
present  work  is  focused  on  further  analysis  into  the  implementation  of  reflection  bound¬ 
ary  conditions.  The  implementation  of  the  boundary  condition  that  was  applied  in  [4] 
is  discussed  further,  and  alternate  higher  order  approach  is  presented.  Cheng  et  al.  [8] 
undertook  a  detailed  study  of  reflection  in  the  level  set  method  to  improve  the  efficiency 
of  the  implementation  proposed  in  [5].  In  that  work,  the  grid  is  augmented  with  the 
boundary  condition,  and  the  CFL  condition  is  adjusted  accordingly  to  maintain  stability. 
In  the  event  that  the  boundary  location  lies  within  Ax  of  a  grid  point,  the  point  clos¬ 
est  to  the  boundary  is  removed.  Flowever,  the  accuracy  analysis  of  the  finite  difference 
WENO  algorithm  relies  on  a  uniform  underlying  grid  (or  at  least  a  smooth  mapping  to 
a  uniform  grid)  [9].  Given  the  errors  that  will  inevitably  be  introduced  into  the  model 
by  uncertainties  in  the  ocean  environment,  it  is  desirable  to  work  with  a  uniform  grid 
in  order  to  maintain  the  convergence  properties  of  WENO  to  the  extent  that  is  possible. 
Another  related  approach  to  reflection  was  proposed  in  [10],  where  a  DG  expansion  was 
applied  to  the  phase  space  variable  and  reflection  boundary  conditions  were  derived  in 
the  coefficient  domain.  The  test  problems  presented  in  this  work  did  not  expose  the  phase 
space  cusp  however,  so  it  is  unclear  how  the  method  would  perform  under  more  general 
initial  conditions.  We  find  that  WENO  methods  perform  very  well  for  such  problems. 
Forrer  and  Jeltsch  [11]  and  also  Forrer  and  Berger  [12]  present  a  boundary  condition  for 
computing  reflections  that  is  very  similar  to  the  higher  order  approach  studied  in  this 
work,  however  they  use  finite  volume  methods,  whereas  the  context  of  the  present  study 
is  WENO  finite  differences.  Finally,  the  work  by  [13]  offers  a  very  thorough  description 
of  the  reinitialization  problem  that  is  applied  here  and  in  [8]  to  obtain  the  signed  distance 
function  to  the  boundary  interface  and  thus  the  boundary  normals. 

1.3  Layout 

The  level  set  equations  and  reflection  boundary  conditions  are  reviewed  in  Section  2, 
and  the  proposed  methods  are  introduced  with  an  example  in  one  spatial  dimension.  A 
discussion  of  stability  and  accuracy  properties  follows  in  Section  3.  Results  for  simple 
test  cases  are  presented  in  Section  4,  and  Section  5  concludes  this  paper. 
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2  Background 

2.1  Level  set  method  for  shallow  water  acoustics 

One  begins  with  Eikonal  equation  for  the  acoustic  phase  function  S(t,x )  as  derived  from 
the  geometric  acoustics  (high  frequency)  approximation  to  the  linear  wave  equation, 

S(f,x)  +  c(x)  |  VS(f,x)|  =0,  (2.1) 

with  x  E  O  C  Rd  and  t  E  R+.  The  speed  of  sound  in  the  water  column  is  given  by  c(x)  > 
0,  and  it  is  assumed  to  be  a  known,  smooth  function.  Consider  the  problem  in  two- 
dimensional  physical  space,  with  x=  (x,z),  and  z  increases  with  water  depth.  The  acoustic 
source  is  located  at  xs  =  (xs,zs).  The  physical  acoustics  model  is  of  an  infinite  line  source 
parallel  to  the  y  axis.  The  same  level  set  model  will  result  from  the  more  commonly 
applied  assumption  of  a  point  source  under  azimuthal  symmetry  under  the  restriction 
xER+,  where  in  this  case  x  is  identified  with  the  range  in  a  cylindrical  coordinate  system. 
Both  geometry  types  are  illustrated  in  Fig.  1.  In  this  work,  we  concentrate  on  the  d  =  2 
case  as  illustrated,  but  the  techniques  can  be  extended  to  higher  dimensional  propagation 
problems. 


Figure  1:  Source  geometry.  The  infinite  line  source  symmetry  is  depicted  on  the  left,  and  on  the  right  the  point 
source  with  azimuthal  symmetry. 

In  the  level  set  method  with  d  =  2,  we  embed  the  source  location  as  the  zero  level  set 
of  a  vector-valued  function  O(f,x,z,0)  :R+  x  Q  x  [—  7i,7r]  — >•  R2.  Here  9  is  the  angle  of  the 
propagation  direction  in  the  reduced  phase  space;  the  magnitude  of  the  generalized  mo¬ 
mentum  vector  is  determined  by  the  sound  speed  c.  It  can  be  shown  that  the  components 
of  the  level  set  function,  <^/,  l  =  1,2,  satisfy  the  following: 

=  £  =  — V-V.  (2.2) 

This  is  the  transport  equation  that  evolves  the  wavefront  in  the  direction  of  its  normal 
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vector  which  is  the  local  ray  direction  given  by  the  velocity  field 

/  ccos(0)  \ 

V=  csin(0)  .  (2.3) 

\Jfsin(e)-Jfcos(fl)/ 

The  objective  is  to  solve  (2.2)  subject  to  initial  conditions,  (pi(0,x,z,9)  =  (p®(x,z,9),  l  =  1,2. 
For  a  point  source  (equivalently  the  projection  of  the  line  source  on  to  the  x—z  plane)  this 
is  simply 

( pi=x-xs ,  cp2=z-zs.  (2.4) 

This  equation  is  to  be  solved  twice  for  two  component  level  set  functions  for  which  the 
zero  level  sets  are  orthogonal  to  one  another  at  the  wavefront.  This  is  consistent  with 
Osher's  formulation  for  geometric  optics  [5]. 

To  compute  numerical  solutions,  apply  the  semi-discrete  approximation  C~CxJrCz+ 
Cq,  where  Ci  f  is  the  product  of  the  relevant  component  of  V  evaluated  on  the  grid,  and  a 
WENO  differential  operator.  The  WENO  operator  may  thus  be  computed  dimension-by- 
dimension.  The  problem  statement  in  (2.2)  and  (2.4)  is  incomplete  without  the  addition 
of  boundary  conditions.  Periodic  boundary  conditions  are  imposed  in  the  9  dimension, 
while  in  x  and  z  the  boundary  conditions  are  determined  by  the  upwind  direction,  dis¬ 
cerned  from  the  sign  of  V,  which  is  fixed  given  a  fixed  propagation  direction  9.  With 
the  upwind  methods,  only  two  types  of  boundary  condition  need  to  be  specified.  An 
inflow  condition  is  needed  to  specify  function  values  flowing  into  the  domain.  For  this, 
it  is  sufficient  to  apply  Neumann  conditions,  for  instance.  Since  the  wavefronts  propa¬ 
gate  away  from  the  source  (and  not  into  the  domain),  small  perturbations  at  inflow  due 
to  assumptions  made  about  function  values  flowing  into  the  domain  will  not  affect  the 
location  of  the  zero  level  sets.  Of  greater  interest  is  the  boundary  condition  applied  along 
a  reflecting  surface  which  occurs  for  example,  at  a  pressure  release  boundary  (air-water 
interface)  or  from  a  hard  rock  surface. 

The  reflection  boundary  condition  at  a  point  along  an  interface  given  by  (x,z&)  is 
derived  from  Snell's  law  and  is  stated  in  [8]  as 


(pit  /  X/Zfo,9yefl )  (p(t/X,Zif/9jnc')  /  (2.5) 

where  9refi  and  9jnc  are  the  reflected  and  incident  angles  respectively.  The  subscripts 
indicating  the  level  set  function  components  have  been  dropped  to  simplify  the  notation. 
The  reflected  angle  satisfies 

@refl=2-@B  @inc  71.  (2-6) 

Elere  0b  is  the  angle  of  the  outward  normal  to  the  reflecting  surface.  When  the  point  (x,Z/,j 
is  in  the  grid  over  the  physical  domain,  applying  the  boundary  condition  is  a  matter 
of  pre-computing  the  incident  angles  and  interpolating  cj)  as  a  function  of  9  at  a  fixed 
point  on  the  boundary  to  approximate  the  value  <p{t,x,Z},,Qinc ).  Note  that  one  must  be 


6 


S.  Martinelli  /  Commun.  Comput.  Phys.,  x  (20xx),  pp.  1-28 


careful  in  implementing  the  interpolation  strategy  as  the  level  set  function  generally  has 
discontinuous  derivatives  at  a  reflection  boundary.  Then  at  each  time  step  we  can  directly 
apply  (2.5)  when  computing  the  spatial  operator.  For  this  purpose,  [4]  presents  a  WENO 
procedure  adapted  to  reconstruction  at  an  arbitrary  point. 

2.2  Computing  boundary  normals 

The  method  proceeds  based  on  the  exposition  in  [8,13,14],  A  level  set  function  represents 
the  boundary  implicitly.  Let  Q  be  the  rectangle  containing  the  portion  of  the  water  col¬ 
umn  being  modeled.  For  convenience,  take  Q  =  { (x,z)  £  [—1,1]  x  [0,2]}.  Overlay  a  grid 


{(XuZj A)}^kl0,  with 


Consider  these  grid  points  to  be  the  centers  of  the  cells  1^= [x,_ 1 /2,-C  i-  i /2] x  [zj-i/irZj+i/i] x 
[#ic— 1/2A+1/2]- 

Assume  the  ocean  surface  is  given  by  z  =  0,  and  the  bottom  is  described  by  z  = 
Zf,(x),  where  Z/,(x)  is  some  smooth  function.  Establish  the  following  notation:  I  = 
{(x,z)  :z  =  z&(x)},  Q+  =  {(x,z)  :z>Zj,(x)},  and  =  O  — T  — Q+.  This  is  illustrated  in 
Fig.  2. 

We  seek  to  compute  0g(x;),  the  angle  of  the  outward  normal  to  the  interface  at  the 
center  of  each  grid  cell  along  the  boundary  7„(^,  z' = 0, 1, ■  ■  •  ,N—  1,  ji,  is  the  index  of  the  grid 
cell  containing  Zg  (x,)  (Fig.  3).  This  computation  is  restricted  to  the  physical  domain  (x-z), 
so  for  now,  ignore  the  k  index. 


2, 
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Figure  2:  Domain  geometry. 
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Figure  3:  A  grid  cell  intersecting  the  boundary. 

A  distance  function,  D(x)  is  defined  by 

D(x)  =  min(|x  — xj|).  (2.7) 

x/cr 

Thus  D(x)  has  the  property  that  D(x,Z;, (x))  =0.  A  signed  distance  function,  d(x)  satisfies 
the  following  additional  properties: 

1.  |d(x)|  =  D(x) 

2.  d(x)  =  0  when  xeT 

3.  d(x)  =  D(x)  when xG Q+ 

4.  d(x)  =  —  D(x)  when  xG  Q 

5.  |  Vrf(x)  |  =  1 

6.  N  =  Vd  is  the  unit  normal  to  the  interface  T. 

In  order  to  obtain  (x;),  we  seek  to  construct  the  signed  distance  function,  d(x,z)  for 
the  interface  z  =  Z\,(x).  Then,  using  Property  6  above,  we  can  compute  N  and  use  the 
result  to  compute  its  angle. 

To  construct  d(x,z),  Sussman,  Smereka,  and  Osher  [15]  suggested  solving  the  follow¬ 
ing  reinitialization  equation 

dT+S(do)  (|  Vrf|  —  1)  =0  (2.8) 

to  steady  state,  where  S(-)  is  a  sign  function  and  do  is  an  initial  value  of  d.  The  term 
reinitialization  equation  comes  from  the  use  of  this  equation  to  reset  the  level  set  function  to 
a  signed  distance  function  in  the  traditional  level  set  method.  Since  the  characteristics  of 
this  equation  propagate  outward  in  the  normal  direction  from  the  interface,  this  method 
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converges  quickly  to  a  signed  distance  function  near  the  interface  (cf.  [13]).  Peng  et  al.  [13] 
refined  the  method  by  replacing  S(do)  in  (2.8)  with  the  smoothed  out  function 


s(d)  = 


d 

yf d2 + d%  A.r4 + d~  Az4 


to  improve  convergence.  We  initialize  with  the  estimate 

d0(x,z)=zb(x)-z. 


(2.9) 


(2.10) 


To  implement  the  method,  we  use  a  fifth  order  WENO  scheme  to  compute  the  spa¬ 
tial  derivatives,  the  same  TVDRK  method  as  is  used  to  solve  the  level  set  equations  to 
discretize  time.  The  partial  differential  equation  (2.8)  is  of  Hamilton-Jacobi  form,  with 
Hamiltonian 

H(\7d,d)  =  s(d)  (|  Vd|  —  1) .  (2.11) 

For  an  upwind  scheme,  we  need  to  define  a  monotone  numerical  Hamiltonian.  We 
choose  Godunov's  method  which  for  this  problem  which  results  in  the  following  dis¬ 
cretization 


where 


H  ( dij  )=Sij(  \]  max((A+)2,(£>-)2)  +max((C+)2,(D-)2)  - 1 

+sjj  ^v/max((A-)2,(B+)2)  +  max((C-)2,(D+)2)-l^, 


A  —  ^  Vx  d^ ,  B  —  ^  Vxdij  / 
C  =  —  Vz  d^ ,  ®  =  ^zdij , 


(2.12) 


V  and  V+  represent  the  one-sided  difference  operators  that  look  to  the  right  and  to  the 
left  respectively,  and 


a+ =  max(a,0),  a  =min(a,0). 

For  an  implementation  with  first  order  time  integration,  the  scheme  then  can  be  written 
down  as 

d:;+1  =  d”-ArH(d”.).  (2.13) 

A  tolerance  may  be  defined  to  determine  when  steady  state  has  been  reached  within 
some  band  of  grid  cells  surrounding  the  boundary,  or  in  practice  one  could  choose  a 
fixed  number  of  iterations.  Central  differencing  may  then  be  used  to  compute  the  value 
of  the  normal  vector  and  extract  0g  for  each  of  the  X{,  i  =  0,-  •  •  ,N— 1. 
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2.3  Reflection 

Fix  (Xi,ek)  and  consider  the  operator  Cz  introduced  in  Section  2.1.  The  WENO  procedure 
yields  reconstructions  <P+(xi,Zj- 1/2 A)  and  (p^+1/2rk~  (p~ {xi,zj+ i/2rdk)  for  j  = 

0,---,N—  1.  The  superscripts  (-)+  and  (•)“  indicate  the  limits  from  the  right  and  left, 
respectively.  This  allows  us  to  write 


Cz  =  Cy  sin0 


<pi,j+l  /2,k  tpij—  1  /2,/c 

h  ' 


where  </>  is  the  numerical  upwind  flux  given  by 


,  / 4>tj+vu'  if  sin(0^  <  °' 

/,/+ 1 /2,/c  otherwise. 


vi,j—\/2,k’ 


(2.14) 


If  Z[,(xj)  =  2,  then  the  grid  point  Z;y-i/2  coincides  with  the  reflecting  boundary  so  (2.5) 
is  implemented  by  setting  (p+N  _i/2k  =  -\/2  N-k- v  f^e  following,  we  consider  the 
situation  where 

Zb(xi)=zjb-1/2+jhr  (2.15) 

with  0  <  7  <  1  for  some  index  fa  €  {0, 1,  •  •  • ,  N — 1 } .  The  one  dimensional  restriction  is  shown 
in  Fig.  4. 


Zb 

Figure  4:  ID  representation  of  the  grid  with  boundary  located  at  Z{,  =  Zy6_i/2  +  7^L 

Two  approaches  are  suggested.  In  the  first,  described  in  Section  2.3.1,  the  location 
of  the  boundary  is  approximated  by  the  nearest  grid  point,  in  Section  2.3.2,  a  different 
approach  is  described  which  assumes  that  the  sound  speed  is  constant  in  a  neighborhood 
of  the  boundary  in  an  effort  to  derive  a  higher  order  approach.  Each  is  described  with 
respect  to  computing  Cz.  The  details  are  similar  for  computing  Cz. 

2.3.1  Method  1:  Boundary  location  approximation 

This  approach  is  straightforward  and  easy  to  implement.  Given  /*,  for  each  xt,  apply 
(2.5)  at  the  location  Zyt_ \/2.  Most  of  the  computation  can  be  performed  prior  to  the  main 
wavefront  solver  loop  and  stored  in  appropriate  data  structures.  The  following  steps  are 
applied  in  the  pre-compute  stage: 
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1.  Solve  the  Hamilton-Jacobi  equation  for  the  signed  distance  function  to  the  boundary  as  described 
in  Section  2.2  and  store  the  values.  Also  store,  for  each  x,,  the  angle  of  the  outward  normal  to 
the  boundary,  0j,(x;). 

2.  Loop  over  each  physical  space  dimension  individually  to  locate  the  indices  4  for  the  left  and  right 
boundaries,  and  jf,  for  surface  and  bottom  boundaries.  Store  the  indices  4  associated  with  each 
j  =  0,  -  ■  ■  ,N  —  1,  and  also  the  indices  4  associated  with  each  i  =  0,-  ■  •  ,N—  1. 

3.  For  each  (4(x;),  loop  over  all  angles  6 j-,  k  =  0,-  ■  ■  ,N  —  1,  to  identify  reflected  angles  (this  may 
be  accomplished  by  testing  the  sign  of  sin(0j-)sin(0B)+cos((4)cos($B))-  Store  the  reflected 
angles. 

4.  For  each  reflected  angle  identified,  compute  and  store  the  associated  incident  angle  using  (2.6). 

2.3.2  Method  2:  Local  ray  approximation 

In  this  approach,  we  combine  the  law  of  reflection  with  the  fact  that  (2.2)  transports  func¬ 
tion  values  along  rays,  which  are  straight  lines  when  the  wave  speed  c  is  constant.  As¬ 
sume  that  c(x,z )  is  a  smooth  function  that  may  be  reasonably  approximated  as  constant 
in  a  neighborhood  of  Zf,(x/),  and  that  Z/,(x)  is  smooth  and  reasonably  approximated  as 
piecewise  linear.  The  time  domain  ray  (characteristic)  equations  for  the  Eikonal  equation 
are 


x(f)  =c(x(t),z(t))cos0(t), 

(2.16a) 

z(t)  =  c(x(t),z(t))sin0(t), 

(2.16b) 

0(f)  =  -^-sin0(f)  —  -j^cos  0(f). 
dx  dx 

(2.16c) 

The  dot  indicates  differentiation  with  respect  to  time  f.  For  small  f,  we  can  expand  0(f)  ~ 
0(0)  =  do,  and  likewise,  using  (2.16), 

x(f)«xo  +  fc(xo,zo)cos0o,  (2.17a) 

z(t)ttz0  +  tc(xo,Zo)sm90.  (2.17b) 

The  assumption  is  that  the  travel  time  along  a  ray  for  which  ||(x,z)  —  (x,Zf,(x))||  <h  is 
small,  so  that  the  rays  can  be  treated  as  straight  lines  in  this  region.  Changing  notation 
for  the  moment,  let  R  =  (x;,z/(i_-i /2),  and  let  P  =  (x/„Z/;)  where  (x;,z/(i _j /2)  are  as  above, 
and  (xf,,Zfo)  is  the  point  of  boundary  interaction  for  the  ray  that  reaches  point  R  at  time  f. 
Let  r  be  the  travel  time  along  that  same  ray  from  the  boundary  interaction  until  it  reaches 
R.  Then  by  property  of  the  transport  equation  for  the  level  set  functions  and  (2.5), 

<p(t,R,0refi)=(p(t-T,Pr0refi)=(p(t-T,P/9inc).  (2.18) 

In  the  absence  of  the  reflecting  boundary,  the  ray  would  continue  to  travel  in  a  straight 
line  to  a  location  R'  £  Q+,  and  by  (2.18),  the  level  set  function  would  satisfy 

<P  ( t,  R,  9  re  fi)=(p(t,R',0inc). 


(2.19) 
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R' 


Zb 


Figure  5:  Law  of  reflection. 


That  is,  if  we  continue  to  propagate  the  level  set  function  into  Q+,  an  exact  boundary 
condition  is  available  when  c  is  (locally)  constant  and  Zf,(x)  is  (locally)  linear  (although 
typically  R'  will  not  lie  in  the  grid  so  interpolation  is  needed).  So  to  implement  the  re¬ 
flection  condition,  extend  the  domain  of  computation  into  Q+  for  directions  Or  incident 
on  the  boundary,  compute  the  locations  R'  for  each  pair  (. Xi,Zjb ),  and  apply  (2.19). 

To  compute  R' ,  consider  the  geometry  illustrated  in  Fig.  5.  According  to  the  law 
of  (specular)  reflection,  ///  =  tjr,  where  ///  is  the  angle  between  the  incident  ray  and  the 
boundary  normal,  and  i]r  is  the  angle  between  the  boundary  normal  and  the  reflected 
ray.  Also,  under  the  constant  wave  speed  assumption,  the  length  of  the  ray  from  P  to 
R'  is  equal  to  the  length  of  the  ray  from  P  to  R,  which  is  cr.  This  implies  that  the  line 
segment  RR'  is  normal  to  the  boundary,  i.e.,  R'  is  just  the  normal  reflection  through  the 
boundary  at  the  point  Q  =  (xi,Zj, (*;)).  Thus,  given  R,  0f-;(x,),  and  assuming  the  boundary 
is  linear  in  the  cell  (*i-1/2,Xi+1/2),  {R')T  =  TRT ,  where 


/  sin2  —  cos2  —  2  cos  sin  0b 

\  —2  cos  0b  sin  0b  —  sin2  0b  +  cos2  0b 


(2.20) 


To  implement  this,  the  following  steps  additional  to  those  listed  in  Section  2.3.1  are 
required: 

5.  Extend  the  domain  into  Q+  to  allow  for  solutions  surrounding  the  image  location.  The  distance 
beyond  the  boundary  of  the  image  location  is  bounded  above  by  yh,  so  the  number  of  grid  points 
to  insert  beyond  the  boundary  can  be  determined  based  on  the  stencil  size  of  the  numerical  solver 
and  interpolation  used,  but  a  minimum  of  2  additional  cells  is  recommended. 


6.  Extend  the  wave  speed  c(x,z)  into  fl+.  If  c  is  constant,  this  is  trivial.  More  generally,  it  is  not 
clear  what  the  best  approach  is.  We  choose  to  set  C;;-  to  be  an  average  of  the  values  obtained 
from  neighboring  points  in  Cl~ ,  fixing  the  value  as  constant  for  (xj,z j)  such  that  the  neighboring 
points  are  all  in  fl+,  and  setting  || \-tj  —  g§ |;-.  =  0.  It  is  less  clear  how  to  extend  c(x,z)  into  Q+ 
if  one  wishes  to  use  a  higher  order  approximation  of  the  ray  in  the  image  space. 


7.  For  each  0j,(x;),  compute  and  store  the  associated  image  mapping  (2.20). 

8.  For  each  0j,(xf),  compute  and  store  the  image  location. 
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2.4  Remark  on  implementation 

The  prototype  code  is  written  in  C  using  the  GNU  Scientific  Library  [16].  In  order  to  store 
all  of  the  precomputed  information,  we  define  the  following  C  structures  in  Listings  1  and 
2. 


Listing  1:  Boundary  gridcell  structure 


typedef  struct  { 

size.t  indx;  II Index  of  this  cell 

int  boundary.indx ;  //Index  of  the  boundary 

double  th_normal;  // Outward  normal  to  the  boundary 

size.t  num_angles;  //#  of  inc/refl  angle  pairs 

size.t  *  refl.indx ;  //Array  of  indices  into  refl.  angles 

size.t  *  non_refl_indx ;  //Complement  of  refl.indx 

gsl_vector  *  inc_angles;  // Associated  incident  angles 

gsl_vector  *  x_im;  //Store  computed  image  locations 

gsl_vector  *  z_im; 

Box  *  neighbors;  //Store  adjacent  points  for  interpolation 
gsl_matrix  *  R;  //Reflection  map 

}  gridcell ; 


Listing  2:  Boundary  Slice  structure 

//  This  is  a  higher  level  structure  to  contain  data  when 
//  there  are  multiple  segments  in  a  vertical  or  horizontal  slice. 

typedef  struct  { 

size.t  Nseg;  //  Index  into  grid  of  this  cell 
gridcell  *  Bleft;  //  An  Nseg— length  array  of  gridcells 
gridcell  *  Bright; 

double  *  zstarl.left;  //  Boundary  limits 
double  *  zstar  1  _right ; 
double  *  zstar2 ; 


}  Slice  ; 

A  Slice  contains  the  boundary  information  for  a  single  dimension  in  the  physical 
space.  Its  members  include  a  gridcell  array  for  each  boundary  in  that  dimension.  If  the 
domain  is  convex,  the  field  Slice .  Nseg,  which  specifies  the  length  of  the  arrays,  equals  1. 
Otherwise,  if  there  are  multiple  segments  (e.g.,  for  an  object  of  interest  inside  the  domain), 
we  store  a  gridcell  for  each  segment,  tied  to  the  relevant  boundary  (left  or  right). 
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3  Analysis 

In  this  section,  we  study  the  proposed  boundary  conditions  more  closely  to  help  quantify 
the  error  that  will  be  imposed  by  the  approximations.  First  we  consider  the  accuracy  of 
the  boundary  conditions  in  terms  of  arrival  time  errors.  Second,  we  look  at  the  stability 
and  show  that  it  is  worthwhile  to  implement  a  higher  order  method  such  as  the  one 
described  in  Section  2.3.2.  It  is  also  worth  noting  that  the  piecewise  linear  assumption  on 
Zf,(x)  (or  equivalently,  piecewise  constant  for  0g(x))  results  in  additional  0(h2)  error. 


3.1  Arrival  time  accuracy 

To  study  the  error  produced  by  this  boundary  condition,  it  is  not  sufficient  to  apply  Taylor 
expansions  since  the  CFL  condition  guarantees  that  (p  is  not  differentiable  in  the  interval 
(zfo  —  h,Z},).  Instead,  we  investigate  specific  errors  that  occur.  In  this  subsection,  we  focus 
on  how  the  wavefront  travel  time  is  affected,  which  is  more  physically  relevant  than 
global  error  in  cp  since  it  provides  the  phase  information  for  the  Eikonal  equation  (2.1). 


3.1.1  Method  1 


Let  t*  be  the  arrival  time  of  the  wavefront  at  the  boundary  Z/:,  and  consider  the  case  where 
c  is  constant,  and  0g(x,)  =  y.  Then  the  incident  angles  lie  in  (0,7r)  and  t*  =  Using 
(2.15),  the  numerical  wavefront  arrival  would  occur  at  time  t*=t*  +  so  that  the  arrival 

time  error,  ef*  is  given  by 


jh 

csinF 


(3.1) 


To  evaluate  further,  consider  0g  =  j.  Then  the  incident  angles  satisfy  0,„c  =  —0ref\.  For 
angles  close  to  normal  incidence,  ef~yr=(D(h)  for  a  first  order  approximation.  The  worst 
case  error  will  occur  for  incident  angles  that  are  close  to  parallel  with  the  boundary.  For 
a  fixed  h,  this  occurs  when  6  =  n(l  —  |).  Then,  ef*  ~  py,  which  appears  to  result  in  an 
order  one  error  in  the  arrival  time  at  these  angles.  In  fact,  we  cannot  guarantee  for  a 
given  refinement  {/7/}^0,  with  hj  — »  0  that  the  corresponding  sequence  {7/}  satisfies 

/— >00 

lim/  ^co7/  =  0.  This  can  be  seen  by  expressing  7 /  in  terms  of  the  first  value  70  as  7 /  = 
if  —  Lif  J  •  The  sequence  {7;}  need  not  converge:  for  example,  take  70  =  7  with  hi  =  2~' . 
On  the  other  hand,  if  the  hj  are  chosen  such  that  there  is  an  index  L  for  which  is  an 
integer  for  l  >=L,  then  ef«  — » 0  as  hj  — >  0.  This  is  essentially  modifying  the  grid  to  include 
the  point  Zf,,  and  that  modification  could  be  different  for  each  fixed  x.  But  in  general,  h— >0 
does  not  imply  ef«  — >  0.  This  error  in  arrival  time  that  is  proportional  to  7  and  dependent 
on  the  arrival  angle  is  actually  a  numerical  consequence  of  the  fact  that  as  h  —r  0,  the 
reflected  angle  approaches  an  angle  that  is  parallel  to  the  boundary  and  is  not  actually 
a  reflected  angle.  Since  this  is  a  consequence  of  the  physics  (i.e.,  the  cusp  in  (p(-,9 ))  we 
would  expect  any  approximation  to  exhibit  this  degradation  at  reflected  angles  that  are 
nearly  parallel  to  Z/:,(.r,). 
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3.1.2  Method  2 

As  mentioned  in  Section  2.3.2,  this  boundary  condition  is  exact  when  c  is  constant  and 
zb(x)  is  linear  in  each  cell  (x;-i/2/*;+l/2),  so  comparable  analysis  to  that  in  Section  2.3.1 
yields  et*  =  0.  Instead,  consider  the  error  made  in  the  approximation  6(t)  ~  do.  The  key 
assumption  is  that  r,  the  travel  time  along  the  ray  that  is  reflected  at  the  boundary  at 
location  P  and  lands  at  location  R  inside  the  domain  (c.f.  Fig.  5),  is  small.  To  investigate 
the  validity  of  that  assumption,  we  solve  for  r  using  the  ray  equations  (2.17)  combined 
with  the  locally  linear  boundary  assumption 

zb  (x ) «  zb  ( Xi ) + cot  (0B  ( Xi ) )  ( x -  Xi ) .  (3.2) 

Proceed  by  taking  P  =  (xb,zb),  where  the  pair  (xb/zb)  satisfy  (3.2),  and  R  =  (xr ,Zr),  where 
|zr  —  zb\  =  7 h  and  0  <  7  <  1.  First  compute  the  travel  time,  Ti,  under  the  approximation 
with  c  locally  constant,  and  compare  to  the  true  travel  time  T2  for  linear  c.  For  T\,  (2.17) 
becomes 


Xr=X(,  +  CTiCOS  9refh 

(3.3a) 

zR=zb+CT1sin9refi, 

(3.3b) 

9(t 1 )  9reji . 

(3.3c) 

Upon  substitution  of  (3.2)  into  (3.3),  we  see  that 

7  h 

c  ( cos  9refj  cot  9b  —  sin  9refi )  ' 

(3.4) 

which  for  9g  =  j  is  exactly  (3.1)  in  magnitude.  Taking  f  f°r  comparison,  the  re¬ 
mainder  term  for  the  small  time  Taylor  expansion  of  0(f)  about  9refi  will  be  proportional 
to 


\Mn)\" 


7// 


sin 


(Orefl)  | 


|isme(f)-|cos«(f) 


(3.5) 


for  some  £€  (0,Ti).  Again,  this  is  order  h  for  reflection  angles  close  to  normal,  and  order 
1  for  angles  near  parallel  with  the  boundary,  with  the  constant  dependent  on  the  sound 
speed  profile,  as  h  — >  0.  Typical  sound  speed  profiles  vary  more  rapidly  in  depth  than  in 
range,  so  we  would  expect  the  [pcos  9  term  to  dominate. 

We  can  now  consider  how  the  approximation  affects  the  travel  time  error  by  looking 
at  the  effect  when  c  is  linear  and  independent  of  x,  since  in  that  case  an  analytical  so¬ 
lution  is  available  so  that  we  may  compute  the  exact  travel  time  T2  and  compare  to  T\ 
to  obtain  an  error  estimate.  Letting  c(z)  =  gz  +  <7  where  g  and  Cq  are  known  constants, 
the  travel  time  along  the  ray  that  leaves  the  boundary  z&  and  terminates  at  Zr  is  given 
by  r=foi  '(HU))dU'  where  s  is  the  arc  length  of  the  ray.  For  the  given  c,  the  travel  time 
is  [17,  p.  210] 


l  +  \/l  —  a2c2 
ac 


c(R) 

c(Zfc) 


T2  = 


(3.6) 
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The  constant  a  is  found  by  applying  Snell's  law  to  the  ray  at  its  turning  point,  yielding 
a  =  CC^I'  ■  Substituting  back  into  (3.6)  and  using  the  fact  that  c(zj,)  =c(R )  +gryh, 


1  /  c(R)+gjh+  sJc(R)2sin2eref,+gjh(2c(R)+g-fh)\ 

Tl  ^l0S\^  c(R)(l  +  |sin0re/;|)  J 

For  reflections  roughly  normal  to  the  boundary,  6refi  m  —  j,  (3.7)  simplifies  to 

1 


t2  = 


=  r1  +  C(yh)2. 


7/7 

gj2h2 

c(R) 

c(R)2 

(3.7) 


(3.8) 


Then,  ef*  =  0(h2)  as  h  — >  0  when  the  approximation  is  applied  for  a  linear  sound  speed 
profile  with  normal  incidence.  Recall  that  in  Section  3.1.1  et*  was  derived  assuming  that 
the  sound  speed  was  constant;  the  error  would  be  larger  had  it  been  assumed  that  c 
was  not  constant  due  to  an  error  in  the  angles.  This  result  suggests  this  method  would 
provide  improved  performance  over  the  location  approximation  method  even  for  the 
variable  coefficient  case. 

Again,  we  would  expect  to  observe  degraded  performance  when  the  angle  the  ray 
makes  with  the  boundary  is  near  parallel.  Given  a  fixed  mesh  spacing  h,  the  worst  case 
would  be  9refi  =  —  n\,  for  which  we  have  ri~  nc[R)  ~0{h2)-  Expanding  (3.7)  about  small 
7 h  leads  to 


r2; 


27 


7Tc(R) 


TC\ 


lc(R) 

27  g 


Vh\- 


7  h 


1  ]  ™(R)\ 

2g7  J 


Then  et*  is  proportional  to 


T 2  — Ti! 


27 


rcc{R) 


27 


nc(R) 


(3.9) 


As  h  — >  0,  et*  — t  —  Ti,  which  is  representative  of  the  fact  that  as  9refi  approaches  parallel 
to  the  boundary,  under  the  linear  c  assumption,  the  ray  has  a  turning  point  at  Z\,  so  that 
it  grazes  the  boundary  but  does  not  reflect.  So  the  arrival  time  error  will  have  the  same 
magnitude  as  that  in  Section  3.1.1  when  a  linear  c  is  approximated  by  constant  c. 


3.2  Stability 

To  simplify  the  analysis,  we  study  the  ID  problem  for  x.  First  we  set  up  the  analysis 
when  the  exact  boundary  condition  is  applied  (i.e.,  the  boundary  location  is  in  the  grid). 
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In  Section  3.3,  we  consider  the  effect  of  errors  incurred  in  approximating  the  boundary 
condition.  In  one  dimension,  there  are  only  two  directions  given  by  6  6  {  —  7r,0} .  Equiva¬ 
lently,  we  solve  the  system 


9U  w  x9U 
aF+A(x)9W 


where 


U  = 


(t,x), 


Initial  and  boundary  conditions  are  required  to  pose  the  problem  completely.  In  this  case, 
we  have 

u(0,x)  =  v(0,x)  =  x  —  x* , 

where  x*  is  a  fixed,  pre-determined  point  in  (0,xg).  It  is  assumed  that  c(x)  >  e  >  0,  so 
that  u(t,x )  represents  a  right-traveling  wave  and  v(t,x)  is  a  left-traveling  wave.  Thus, 
boundary  conditions  need  to  be  imposed  at  x  =  0  for  u  and  at  x  =  Xg  for  v.  Assume  for 
now  that  reflection  only  occurs  at  x  =  xB,  and  apply  an  inflow  condition  at  x  =  0: 


6  ;)?«-• 

(_°1  °)u(f,xB)  =  0. 


The  two  equations  are  now  coupled  by  the  reflection  boundary  condition. 

To  further  simplify  the  analysis,  consider  the  constant  coefficient  problem  in  the  half 
plane  x£(0,oo).  We  will  use  standard  notation  for  k=At,  A=|,  and  E-i  is  the  shift  operator: 
Ei(U(x))  =  U(x+jh).  Apply  the  first  order  scheme  with  reflection  at  x  =  0: 


U'!+1  = 


cA  0\  i 


0  0 


E~i  +  (l-cA)I  + 


0  0 
0  cA 


U/ 


=  QU”=(a_1E-1+A0E°+A1E1)u”,  j  =  1,2,--, 


u'!+1  — 
u0  — 


1-cA 
1  — cA 


0  cA 
0  cA 


TI'! 

u0. 


(3.10) 

(3.11) 


We  first  study  the  effect  of  only  the  boundary  condition  on  the  left,  using  the  techniques 
presented  in  [18]  and  [19]. 

Under  a  few  assumptions  on  Q  (which  are  satisfied  for  this  scheme;  refer  to  [18]),  Q 
has  two  linearly  independent  right-going  and  two  linearly  independent  left-going  solu¬ 
tions  when  |£|  >  1.  It  is  established  in  [18]  that  the  number  of  independent  right-going 
solutions  is  determined  by  the  size  of  the  system  multiplied  by  the  number  of  stencil 
points  in  Q  to  the  left  of  center,  and  likewise  for  the  left-going  solutions  with  respect  to 
the  number  of  stencil  points  to  the  right  of  center.  We  now  consider  solutions  in  the  form 

UJ  =  Uk’C,  K,£eC-{0},  (3.12) 
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with  k  =  e‘th  and  £  =  e*  .  The  dispersion  relation  relates  go  to  £  for  waves.  Analogously, 
we  follow  [19]  and  consider  the  numerical  dispersion  relation  for  £  in  terms  of  k.  This  is 
given  by  substituting  (3.12)  into  Q  and  solving.  We  find  that  the  resulting  system  satisfies 


£-t-i+‘A  0 

0  [cA(k  —  1)  +  1], 


The  matrix  on  the  left  hand  side  is  singular  when 

C  =  cA— -  +  1,  cA(k  — 1)  +  1.  (3.13) 

The  first  solution  is  for  the  right-going  solution  component  u(t,x)  and  the  second  is  for 
the  left-going  v(t,x).  In  the  subsequent,  we  will  use  }i  instead  of  k  for  the  left-going 
solutions  and  retain  k  only  for  the  right-going  solutions,  so  the  second  dispersion  relation 
in  (3.13)  becomes 

£  =  cA(p  — 1)  +  1. 

We  now  seek  solutions  to  the  scheme  (3.10)  that  also  satisfy  (3.11)  and  have  the  form 


2  2 

Uf  =  C  E  UmPmK’m+C  E  / (3.14) 

m= 1  m=l 


Continuing  to  follow  Trefethen  [19],  apply  the  boundary  condition  (3.11)-(3.14)  and  col¬ 
lect  terms  in  a  and  /3  to  arrive  at  an  expression  of  the  form 

E<o(:;)=D©(ft). 

Our  concern  is  solutions  for  which  the  coefficients  a  cannot  be  fully  specified  in  terms  of 
/3  when  |£|  >  1,  i.e.,  |£|  >  1  for  which  E(£)  is  singular.  We  see  that 

Fm=/W  +  [cA(l-Kl)-l]^  C$2  +  [cA(l  — k2)  —  1]  p|\ 

V  (£+cA(l-Kl)-l)P2  (£+cA(1-k2)-1)P2  )' 

/-CKl  +  [cA(/q-l)  +  l]R2  -CRi  +  [cA(^2-l)  +  l]R2\ 

{  (l-C+cAfa-l  ))&i  (1-^+cA(F2-1))R22  )' 

so  E(£)  is  singular  when 

£det(P)(£+cA(l-Jc)-l)=0.  (3.15) 

Here,  P=  [PiP2]  and  we  have  used  the  fact  that  the  dispersion  relation  for  the  scheme 
produces  only  one  solution  for  right  going  waves  so  that  K\  =  k2  =  k.  In  other  words, 
instability  is  present  (for  £  G  £—  {0})  when  P  is  singular,  implying  that  there  are  not  two 
linearly  independent  right  going  waves,  or  £= 1  —  cA(l  —  k) .  In  the  case  that  P  is  singular, 
the  problem  is  one  dimensional  and  the  condition  is  trivially  satisfied  using  the  conse¬ 
quence  of  the  dispersion  relation  that  jc=  A,  with  a  reflection  coefficient  E(^)“1D(^)  =  1. 
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On  the  other  hand,  if  P  has  full  rank,  then  £(£)  is  singular  only  when  £=1  — cA(l— k).  Ap¬ 
plying  once  again  the  dispersion  relation,  there  are  two  possibilities,  £  =  1  or  £  =1  —  2cA. 
If  £  =  1,  then  k  =  j-t  =  1  yielding  the  solution  U(f,x)  =  constant.  So  the  approximation  is 
stable  as  long  as  cA  <  1  since  then  1 1  —  2c A  |  <  1. 

Again  as  in  [19],  we  can  define  the  reflection  coefficient  matrix  Ar(£)  =  E_1(£)D(£) 
when  E  is  nonsingular.  It  is  straightforward  to  compute 


which  is  independent  of  £.  This  matrix  is  actually  singular  and  yields  solutions  a 2  cxfti 
which  suggests  that  the  analysis  for  the  case  P  is  singular  is  the  correct  analysis,  with 
reduced  reflection  coefficient  of  one. 

In  a  similar  fashion,  it  can  be  shown  that  the  problem  for  a  reflection  boundary  at  x=l 
with  boundary  condition 


U 


n+1 _ 

N  ~ 


fl-c  A 
yl  — cA 


U 


n 

N 


is  also  stable  with  one  linearly  independent  solution  in  each  direction,  and  reflection 
coefficient 

Ar(0  =  k2N. 

This  is  stable  since  for  |£|  >  1,  |k|  =  |  ^-i+cA  I  —  1-  By  Proposition  6  in  [19],  the  model  is 
stable  in  the  sense  defined  in  [18]  because  the  interfaces  at  x  =  0  and  at  x  =  1  are  both 
individually  stable. 


3.3  Effect  of  perturbations 

The  proposed  boundary  conditions  rely  on  an  approximation  thus  an  error  is  imposed  at 
each  reflecting  boundary  at  each  time  step.  In  this  section,  we  study  how  the  error  made 
at  the  boundary  propagates  in  the  one  dimensional  case  above.  Consider  (3.10)  with 
initial  condition  Uy  =0  and  perturbed  boundary  condition  Uq=Vq+6,  where  e  is  the  0(h ) 
or  0(h 2)  approximation  error  for  the  boundary  condition.  We  also  have  the  reflection 
boundary  condition  at  x  =  1  for  v'  VN  =  UN ■  We  consider  the  u'j  and  v'j  components 
separately  in  the  scheme  (3.10)  to  construct  the  discrete  error  equation. 

We  seek  to  solve  the  following  difference  equation: 

u'j  =  (1  —  cA)u"-1  +c\uJZi  (3.16) 

subject  to  the  initial  condition  iy°  =  0  and  boundary  condition  /y'j  =  e.  The  transport  prop¬ 
erty  further  implies  that 


u'l  =  0  for  j  >  n. 


(3.17) 
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We  use  generating  functions  to  solve  the  difference  equation.  First,  let 


OO  11  —  1 

Un(x)  =  YjU'jX ’  =  Yj  U'jxi' 
j= 0  7=0 


(3.18) 


where  the  second  equality  is  due  to  (3.17).  Now  multiply  (3.16)  by  x>  and  sum  over  j  to 
find 


n—2 


n  vk+l 


u„+1(x)  =  (l-cA)un(x)+  Y  cAulx‘ 

k=-l 

=  (1  —  cA)un  (x)  +cAxun  (x)  +cAun_  1 
=  (1  — cA+cAx)m„(x). 

y" 


Now  let  u(x,y )  =  Y/n=oun (*)  n\<  multiply  (3.19)  by  and  sum  over  n: 

oo  1/H—l  oo  n 

Yu»(x )  =  Y(1~cA+cAx)u"(x)ir- 

n= 1  {  1h  n= 0  n- 

Or, 


9 u(x,y ) 
dxJ 


■  (1  — cA  +  cAx)i<(x,y), 


which  has  the  solution  u(x,y )  =  C(x)^1  cX+cXx)v .  To  find  C(x), 

(X)  -1 

w(x,1)  =  C(x)£’1-cA+cAx-^  '' 

■un(x)  =  C(x)  (1  — cA+cAx)". 


'=Yu»(x) 

n= 0 
\n 


Taking  n  =  1, 
Thus, 

and 


C(x)  (1  — cA+cAx)  =  !<q  +  !/|x T- u^x^  —  =  e 


w(x,y)  =  - 


—  cA  +  cAx 


p(l— cA+cAx)y 
^  / 


u„(x)  =e(l  —  cA+cAx) 


n — 1 


(3.19) 


(3.20) 


(3.21) 


Comparing  the  above  expression  to  (3.18)  and  applying  the  binomial  theorem,  we  see 
that 

' £(n^1)  (1  — cA)n_,_1  (cA);,  ;  =  0,l,---,n  — 1, 

0,  otherwise. 


«"  =  ■ 


(3.22) 


Next,  solve  the  difference  equation  produced  by  (3.10)  for  vj  subject  to  a  perturbation  at 
the  boundary  using  the  same  techniques.  The  difference  equation  is 

v^_j  =  {l-cA)v^-}j+cAvnN}j+1, 


(3.23) 
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subject  to  initial  condition  v%_j  =  0,  and  boundary  condition  =  e,  n  =  1,2,  ■■  ■  It  is  a 
consequence  of  the  transport  equation  that  v^_  -=0  for  j>n,  and  we  can  use  this  property 
in  addition  to  the  boundary  condition  to  derive  the  useful  identity 


We  use  the  expansion 


to  find 


^_„+1  =  e(cA)”-1. 

(3.24) 

n— 1 

Mx)='EVN-jX 

;=o 

(3.25) 

*1 _ )  cA 

N+l  (x)  =  (1  -  cA)vn  (x)  +  —  £  v"  :x~’ +cAv'^+1  -  —  V^_n+1X 

X  j= 0 


-(«-!) 


Using  (3.24)  and  that  =  0  yields  the  difference  equation  in  n: 

Vn+ l(x)  =  ^1  — CA+  Vn(x)-£  ■ 

Again,  let  v(x,y)  =  J^=0vn(x)lj^.  Then  the  resulting  differential  equation  is 


dv(x,y) 

dlJ 


cA\ 


(  cA 


1  — cA+—  )  v(x,y)—eex p  —  y 


x  j 


V  x 


(3.26) 


Eq.  (3.26)  produces  the  homogeneous  solution 

V\{x,y)  =  C(x)exp 


_  ,  cA 

\-CA  +  —)y 


and  particular  solution  V2(x,y)  =  D(x)exp  [^y] .  Upon  substitution  into  (3.26)  we  find 
that  D(x)  =  i^x,  so  that  the  generating  function  for  the  coefficients  v"N  ■  is 


Now 


v(x,y)  =  C(x)exp 


v(x,l)  =  C(x)exp 


.  .  cA 

i-cA+t  jy 


_  ,  cA  . 

l-cA+-  )y 


+ 


1-cA 


exp 


cA 


+ 


1-cA 


exp 


-y 


cA 

— y 

X 


n= O' 


cA\ 1 


=  E^i  C(x)  h-cA+-  +T 

„nn\  V  V  x  )  1 


cAV 
—  cA\x  ) 


CO  -1 

=  LMX)^ 

n= 0 
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I’”w=cw(i-cA+f)  +t7a(t)  ' 


Considering  n  =  1, 


Mx)  =  C(x)  (1-cA  +  y)+  I 

n — 1 

=  'EVN-iX~1  =  VN  =  £' 


£  cA 
-d  i 


which  leads  to  the  result 


C(x)=£- 


1-cA-f 


(l-cA)(l-cA+f)' 

Substituting  this  back  into  the  expression  for  vn  (x),  we  have 

!’"W  =  IVrt  (1-cA-f)(1-cA+f)  +(v)  '  (3'27) 

Finally,  in  order  to  compare  (3.27)  to  the  coefficient  terms  in  (3.25),  expand  using  the 
binomial  theorem: 


Vn(x)=£  (  1  — cA  + 


x(l  — cA)  V  / 


(1  — cA) 


n—\—j  ( 


— (-V 

1  — cA  1  X  / 


—  €  [  1  —  cA  + 


n- 1  j!-l 


n—l—j  f  A  \  ^ 


(l-cA)f 


cA\° 

“)  +,5 


n— 1\  In  — 1 


(i— cA  r'-i(c4)' 


=  Ee(l-cA)’ 

;=o 


n  —  2j\  f  n\  f  cA 


That  is. 


4_;  = 


»  /  V/7  Vi-cA 


e(l -cA)n_1  (^ir^) (") ( Aib  j=0,l,-,»-l, 
0,  otherwise. 


(3.28) 


We  can  now  apply  (3.22)  and  (3.28)  to  explore  how  the  error  produced  at  the  left  boundary 
propagates.  The  perturbation  at  x=0  reaches  the  right  boundary,  x=l,  at  time  t=(N+l)k, 
with  value  u ^  !  1  =  (cA) N  e.  The  reflection  boundary  condition  at  x  =  1  is  u^+1  =  !/^+1  = 
(cA)n£.  Assuming  the  boundary  condition  on  the  right  is  exact,  then  the  perturbation 
returns  to  x  =  0  as 


N+ 1 


)  (1-cA)" 


(3.29) 


22 


S.  Martinelli  /  Commun.  Comput.  Phys.,  x  (20xx),  pp.  1-28 


Upon  the  next  reflection  at  x=0,  this  means  there  is  an  additive  error,  u^N+1  =€+VqN+1  .  It 
is  evident  that  the  error  continues  to  propagate  through  the  domain.  We  can  use  Stirling's 
approximation 

, _  /  7U\  N 

N'.^VlrcNij)  (3.30) 

to  get  a  better  sense  of  the  size  of  the  error  for  large  n,  with  n  expressed  as  multiples  of 
N,  so  this  is  also  a  mesh  refinement  study.  Substituting  (3.30)  in  (3.29), 


v2N+1 

v0 


4(cA)2(1-cA) 


lN 


2 

(n+i)Vtzn' 


(3.31) 


The  limiting  behavior  of  (3.31)  is  overwhelmed  by  the  exponential  term  in  the  numerator. 
Write  g(c\)  =  (cA)2(l  —  cA).  The  base  of  the  exponential  term  is  then  4g(cA).  Since  cA<l 
for  stability,  we  have  g(c\)  Thus  4g(cA)  <  1,  and  although  the  error  persists,  it 

decreases  exponentially  as  h  — >  0  so  it  does  not  have  a  significant  effect  on  either  method's 
convergence  rate.  To  be  thorough,  we  computed  the  error  multiplier  with  cA  =  |  (this 
value  of  cA  maximizes  g{cX)  over  (0,1))  for  N  =  10,  N  =  20,  and  N  =  40.  The  added  errors 
were,  respectively,  e'8.554 x  10  5,  eT.7026xl0  7 ,  and  el.76xl0“12,  where  the  result  for 
N  =  40  is  based  on  an  approximation  for  large  N. 


3.4  Remark 

In  Sections  3.2  and  3.3,  the  analysis  is  limited  to  the  one-dimensional  case.  The  fact  that 
in  this  case  the  reduced  phase  space  is  discrete,  consisting  of  only  two  directions,  permits 
this  type  of  study.  Local  stability  in  two  dimensions  follows  for  any  pair  of  incident 
and  reflected  angles  on  the  surface  by  coordinate  transform.  An  extension  to  fully  two- 
dimensional  propagation  is  left  to  future  work,  although  the  experimental  results  for  the 
two-dimensional  propagation  problem  suggest  the  method  is  stable. 


4  Results 

4.1  Validation  of  arrival  time  analysis 

Computational  results  for  the  arrival  time  errors  are  shown  in  Figs.  6  and  7  below.  In  both 
cases,  the  scenario  is  for  an  upward  sloping  bottom  with  =  -j.  Fig.  6  is  for  a  constant 
wave  speed  c  =  Co  =  1.5  km/sec,  and  for  Fig.  7,  the  wave  speed  is  given  as  c(z)  =0.5z+Co. 
In  each  figure,  the  subplots  on  the  left  display  the  error  computed  at  time  T  =  0.3  seconds 
along  the  wavefront  as  parameterized  by  local  (reflected)  ray  direction  6.  The  subplots 
on  the  right  hand  side  show  the  error  as  a  function  of  grid  size  N  at  the  angle  0  =  — 
which  is  normal  to  the  sloping  bottom. 

From  Fig.  6,  it  appears  that  the  convergence  rate  for  Method  2  is  greater  than  for 
Method  1.  Although  we  claim  that  Method  2  is  exact  for  this  scenario,  we  would  not 
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Method  1  Error  at  0  =  -2ju/3 


Method  2  Error  at  0  =  -2ju/3 


Figure  6:  Arrival  time  errors  with  c  =  1.5  km/sec. 


Figure  7:  Arrival  time  errors  with  c  =  0.5z  +  1.5  km/sec. 


expect  to  observe  convergence  that  is  any  faster  than  0(h2)  due  to  coupling  with  the 
errors  in  the  PDE  solver  and  wavefront  extraction  routine.  We  see  that  the  convergence 
rate  for  Method  1  is  roughly  0(h)  and  for  Method  2,  about  0(h2).  The  error  levels  off 
for  Method  2  at  about  10  5,  this  is  likely  due  to  time  discretization  error.  In  fact,  for 
At  =  0.0087  with  N  =  160,  there  is  an  error  on  the  order  (At)  due  to  the  third-order 
Runge-Kutta  time  integration.  This  leads  to  an  error  (At)  x  Nf  =  2.2x  10  5  where  Nt  =  34 
is  the  number  of  times  steps  processed  to  compute  the  arrival  time  error. 
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Fig.  7  confirms  that  for  the  linear  sound  speed,  the  convergence  rate  slows  as  the 
reflected  angle  becomes  closer  to  parallel  with  the  bottom  in  both  Method  1  and  Method 
2.  The  rate  close  to  9  =  —  ^  is  roughly  0(h)  for  Method  1.  It  is  less  clear  for  Method 
2,  but  the  results  appear  to  be  converging  faster  than  0(h )  at  the  same  angle.  For  both 
plots  the  errors  noticeably  dip  near  0  =  —  j.  At  this  angle,  the  velocity  component  in  the 
x-di  recti  on  is  zero,  reducing  the  dimensionality  of  the  problem. 

4.2  Wavefront  examples 

To  provide  some  visual  results,  we  present  a  few  examples  of  wavefronts  extracted  from 
the  level  set  functions  in  the  presence  of  non-conforming  grid  geometry  with  reflections. 
The  visualizations  in  this  section  were  produced  using  the  Mayavi  program  [20]. 

4.2.1  Comparison  between  Method  1  and  Method  2 

Fig.  8  shows  the  results  from  a  40x40x40  grid  over  Q  =  {(x,z)  E  [ — 0.5,0.5]  x  [0,1]}  and 
a  sloping  bottom  having  6b  =  j  radians.  A  constant  sound  speed  c  =  1.5  km/ sec  is  used 
here.  At  this  resolution,  it  is  clear  that  Method  2  is  preferred,  as  the  artifacts  of  the  bound¬ 
ary  representation  used  in  Method  1  are  apparent.  Note  that  although  the  representa¬ 
tion  of  the  bottom  location  is  piecewise  constant,  the  proper  normals  are  applied.  When 
the  grid  size  is  increased  to  160  x  160  x  160,  the  wavefront  produced  using  the  Method  1 
boundary  condition  is  converging,  though  very  slowly  at  angles  close  to  parallel  with  the 
boundary. 


t  =  0.30  sec 


t  =  0.30  sec 


(a)  Method  1  (b)  Method  2 

Figure  8:  Wavefront  comparison  with  N  =  40  at  0.3  seconds,  c  =  1.5  km/s,  = 


4.2.2  Examples  with  non-linear  geometry 

Figs.  10  and  11  are  examples  of  the  types  of  geometry  for  which  Method  2  provides  nice 
results.  Again,  a  constant  sound  speed  c  =  1.5  km/ sec  is  used  for  simplicity.  The  domain 
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t  =  0.30  sec  t  =  0.30  sec 


(a)  Method  1  (b)  Method  2 

Figure  9:  Wavefront  comparison  with  N  =  160  at  0.3  seconds,  c  =  1.5  km/s,  6g  =  ^. 

for  each  is  Q  =  { (x,z)  £  [—1,1]  x  [0,1]}.  Each  subfigure  is  a  snapshot  of  the  wavefront 
at  a  few  times,  t  =  0.2,0.4,0. 6,0.8.  In  Fig.  10,  a  modified  Gaussian  function  is  applied  to 
represent  an  undersea  canyon,  and  in  Fig.  11,  a  seamount  is  represented  by  a  Gaussian 
function.  Note  that  the  Gaussian  function  is  not  a  limitation,  just  a  convenient  choice. 


-I  OO  OOO  100 

x  (km) 


(a) 

t  =  0.60  sec 


t  =  0.40  sec 


(b) 


t  =  0.80  sec 


(c) 


(d) 


Figure  10:  Wavefront  snapshots  with  z^(x)  =  l 
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t  =  0.20  sec 


t  =  0.40  sec 


(C)  (d) 

Figure  11:  Wavefront  snapshots  with  Zf,(x)  =  1  — 0.3exp  [— 12.5x2] . 

Due  to  the  multiple  reflections  produced  by  the  canyon  example,  we  were  not  able  to 
produce  reasonable  results  using  Method  1. 

5  Conclusion 

Level  set  methods  provide  an  alternative  framework  for  numerical  solutions  to  the  high 
frequency  approximation  to  the  wave  equation.  A  fixed-grid  point  of  view  is  advanta¬ 
geous  for  medium  range  (about  1  kilometer),  high  frequency  modeling.  In  particular,  it 
is  very  difficult  to  compute  accurate  eigenray  solutions  for  active  sonar  in  shallow  water 
environments  with  range-dependent  bathymetric  features.  Reflections  are  an  important 
consideration  for  shallow  water  acoustics,  and  a  reasonable  approximation  at  high  fre¬ 
quencies  where  the  amount  of  energy  transmitted  into  the  sea  bottom  is  negligible. 

We  have  presented  a  modified  level  set  method  directed  toward  computational  un¬ 
derwater  acoustics  modeling.  Two  options  are  proposed  for  boundary  conditions,  the 
first  (Method  1)  using  approximate  boundary  locations  to  maintain  a  uniform  grid,  and 
the  second  (Method  2)  taking  advantage  of  the  transport  property  of  the  level  set  equa¬ 
tion  and  the  law  of  reflection.  Although  both  yield  correct  arrival  times  in  the  large  N 
limit  at  most  angles,  the  method  based  on  the  law  of  reflection  is  shown  to  converge 
about  one  order  faster  than  the  method  using  an  approximate  location.  An  analysis  of 
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the  effect  of  the  boundary  condition  on  arrival  times  is  supported  by  the  results  in  Section 
4.1.  The  computational  examples  in  Figs.  10  and  11  above  show  that  we  are  able  to  com¬ 
pute  very  nice  wavefronts  on  less  restrictive  geometry  types  produced  inadequate  results 
using  Method  2.  It  is  straightforward  using  the  grid  structures  presented  to  extend  the 
approach  to  any  smooth  boundary  representation,  including  bathymetry  data  that  have 
been  smoothed  and  interpolated  to  yield  a  functional  representation.  Thus,  the  tools  are 
available  to  study  more  practical  propagation  problems  using  the  level  set  method. 
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